{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#根据样地点提取纹理特征\n",
    "from osgeo import gdal\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import os\n",
    "import osr\n",
    "from pandas import set_option"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class change_coordinate():\n",
    "    def __init__(self, dataset):\n",
    "        self.dataset = dataset\n",
    "\n",
    "    def getSRSPair(self,dataset):\n",
    "        '''\n",
    "        获得给定数据的投影参考系和地理参考系\n",
    "        :param dataset: GDAL地理数据\n",
    "        :return: 投影参考系和地理参考系\n",
    "        '''\n",
    "        prosrs = osr.SpatialReference()\n",
    "        prosrs.ImportFromWkt(self.dataset.GetProjection())\n",
    "        geosrs = prosrs.CloneGeogCS()\n",
    "        return prosrs, geosrs\n",
    "\n",
    "\n",
    "    def lonlat2geo(self, lon, lat):\n",
    "        '''\n",
    "        将经纬度坐标转为投影坐标（具体的投影坐标系由给定数据确定）\n",
    "        :param dataset: GDAL地理数据\n",
    "        :param lon: 地理坐标lon经度\n",
    "        :param lat: 地理坐标lat纬度\n",
    "        :return: 经纬度坐标(lon, lat)对应的投影坐标\n",
    "        '''\n",
    "        prosrs, geosrs = self.getSRSPair(self.dataset)\n",
    "        ct = osr.CoordinateTransformation(geosrs, prosrs)\n",
    "        coords = ct.TransformPoint(lon, lat)\n",
    "        return coords[:2]\n",
    "\n",
    "\n",
    "    def geo2imagexy(self, x, y):\n",
    "        '''\n",
    "        根据GDAL的六 参数模型将给定的投影或地理坐标转为影像图上坐标（行列号）\n",
    "        :param dataset: GDAL地理数据\n",
    "        :param x: 投影或地理坐标x\n",
    "        :param y: 投影或地理坐标y\n",
    "        :return: 影坐标或地理坐标(x, y)对应的影像图上行列号(row, col)\n",
    "        '''\n",
    "        trans = self.dataset.GetGeoTransform()\n",
    "        a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])\n",
    "        b = np.array([x - trans[0], y - trans[3]])\n",
    "        return np.linalg.solve(a, b)  # 使用numpy的linalg.solve进行二元一次方程的求解\n",
    "\n",
    "    def lonlat2rowcol(self,lon,lat):\n",
    "        '''\n",
    "        根据经纬度转行列公式直接转换为行列\n",
    "        '''\n",
    "#         tp = self.lonlat2geo(lon,lat)\n",
    "        geo = self.dataset.GetGeoTransform()\n",
    "#         row = int((tp[0] -geo[0]) / geo[1]+0.5)\n",
    "#         col = int((tp[1] - geo[3]) /geo[5]+0.5)\n",
    "        row = int((lon -geo[0]) / geo[1]+0.5)\n",
    "        col = int((lat - geo[3]) /geo[5]+0.5)\n",
    "        \n",
    "        return row,col"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "class define_window():\n",
    "    '''\n",
    "    :param w 定义窗口大小\n",
    "    :param center_row 中心点行号\n",
    "    :param center_col 中心点列号\n",
    "    '''\n",
    "    def __init__(self,w):\n",
    "        self.w = w\n",
    "    def window_upleft_rowcol(self,center_row,center_col):\n",
    "        upleft_row = center_row - (self.w-1)/2\n",
    "        upleft_col = center_col - (self.w-1)/2\n",
    "        return upleft_row,upleft_col"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class make_feature_names():\n",
    "    '''\n",
    "    根据波段编写特征名称，返回特征名称列表\n",
    "    '''\n",
    "    def __init__(self,dataset):\n",
    "        self.nb = dataset.RasterCount\n",
    "    def feature(self,feature_list):\n",
    "        names = []\n",
    "        for i in range(self.nb):\n",
    "            for j in feature_list:\n",
    "                names.append('{}{}{}'.format(j,'_',i))\n",
    "        return names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "if __name__ == '__main__':\n",
    "    '''\n",
    "    把图像与坐标放到一个文件夹下\n",
    "    '''\n",
    "    img_dir = r'./使用数据'\n",
    "    out_path=r'./输出数据'\n",
    "    gdal.AllRegister()\n",
    "    img = gdal.Open(os.path.join(img_dir,'500_b0_win7_texture.tif'))\n",
    "    ds = pd.read_excel(os.path.join(img_dir,'point.xls'))\n",
    "    ns = img.RasterXSize\n",
    "    nl = img.RasterYSize\n",
    "    run_change_coordinate = change_coordinate(img)#调用坐标转换函数\n",
    "    w = 7 #窗口大小\n",
    "    run_define_window = define_window(w)#调用窗口定义函数\n",
    "\n",
    "    run_make_feature_names = make_feature_names(img)#调用特征名称函数\n",
    "    names = [ 'mean_1','variance_1','homogeneity_1','contrast_1','dissimilarity_1','entropy_1','sencond_moment_1','correlation_1',\n",
    "              ] \n",
    "\n",
    "    lon,lat = ds.iloc[:,1].values,ds.iloc[:,2].values\n",
    "    '''\n",
    "    定义输出列表\n",
    "    :all_out输出每个窗口下所有特征的值\n",
    "    :all_mean输出每个窗口下所有特征的平均值\n",
    "    :all_std输出每个窗口下所有特征的标准差\n",
    "    '''\n",
    "    \n",
    "    all_out = []\n",
    "    all_mean = []\n",
    "    all_std = []\n",
    "    for i in range(len(lon)):\n",
    "        ilon,ilat = lon[i],lat[i]\n",
    "        ix,iy = run_change_coordinate.lonlat2rowcol(ilon,ilat)\n",
    "        if ix<0 or ix >ns-1 or iy <0 or iy >nl-1:\n",
    "            print('not in the image: '+str(ds.iat[i,0].value))\n",
    "        upleft_x,upleft_y = run_define_window.window_upleft_rowcol(ix,iy)\n",
    "        ref = img.ReadAsArray(int(upleft_x),int(upleft_y),w,w)\n",
    "\n",
    "        if len(ref.shape) == 3:\n",
    "            df = np.zeros((w*w,len(names)))         \n",
    "            for j in range(len(names)):\n",
    "#                 print(j)\n",
    "                df[:,j] = list(ref[j].flatten())\n",
    "            df = pd.DataFrame(df,columns=names)\n",
    "        else:\n",
    "            df = pd.DataFrame(ref.flatten())\n",
    "        description = df.describe()\n",
    "        df_mean = description.iloc[1,:]\n",
    "        df_std = description.iloc[2,:]\n",
    "\n",
    "        all_out.append(df)\n",
    "        all_mean.append(df_mean)\n",
    "        all_std.append(df_std)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "out = pd.concat(all_out)\n",
    "out_mean = pd.concat(all_mean)\n",
    "out_std = pd.concat(all_std)\n",
    "out.to_csv(os.path.join(out_path,'out.csv'))\n",
    "out_mean.to_csv(os.path.join(out_path,'out_mean.csv'))\n",
    "out_std.to_csv(os.path.join(out_path,'out_std.csv'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>mean_1</th>\n",
       "      <th>variance_1</th>\n",
       "      <th>homogeneity_1</th>\n",
       "      <th>contrast_1</th>\n",
       "      <th>dissimilarity_1</th>\n",
       "      <th>entropy_1</th>\n",
       "      <th>sencond_moment_1</th>\n",
       "      <th>correlation_1</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>146.805557</td>\n",
       "      <td>0.679847</td>\n",
       "      <td>1.840844</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.500000</td>\n",
       "      <td>0.750000</td>\n",
       "      <td>0.459098</td>\n",
       "      <td>0.172068</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>146.902771</td>\n",
       "      <td>0.690405</td>\n",
       "      <td>1.839048</td>\n",
       "      <td>0.416667</td>\n",
       "      <td>0.416667</td>\n",
       "      <td>0.791667</td>\n",
       "      <td>0.562930</td>\n",
       "      <td>0.178627</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>146.944443</td>\n",
       "      <td>0.684935</td>\n",
       "      <td>1.815372</td>\n",
       "      <td>0.388889</td>\n",
       "      <td>0.388889</td>\n",
       "      <td>0.805556</td>\n",
       "      <td>0.585526</td>\n",
       "      <td>0.189043</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>146.944443</td>\n",
       "      <td>0.684935</td>\n",
       "      <td>1.815372</td>\n",
       "      <td>0.388889</td>\n",
       "      <td>0.388889</td>\n",
       "      <td>0.805556</td>\n",
       "      <td>0.585526</td>\n",
       "      <td>0.189043</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>147.000000</td>\n",
       "      <td>0.645497</td>\n",
       "      <td>1.621216</td>\n",
       "      <td>0.277778</td>\n",
       "      <td>0.277778</td>\n",
       "      <td>0.861111</td>\n",
       "      <td>0.666667</td>\n",
       "      <td>0.260031</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       mean_1  variance_1  homogeneity_1  contrast_1  dissimilarity_1  \\\n",
       "0  146.805557    0.679847       1.840844    0.500000         0.500000   \n",
       "1  146.902771    0.690405       1.839048    0.416667         0.416667   \n",
       "2  146.944443    0.684935       1.815372    0.388889         0.388889   \n",
       "3  146.944443    0.684935       1.815372    0.388889         0.388889   \n",
       "4  147.000000    0.645497       1.621216    0.277778         0.277778   \n",
       "\n",
       "   entropy_1  sencond_moment_1  correlation_1  \n",
       "0   0.750000          0.459098       0.172068  \n",
       "1   0.791667          0.562930       0.178627  \n",
       "2   0.805556          0.585526       0.189043  \n",
       "3   0.805556          0.585526       0.189043  \n",
       "4   0.861111          0.666667       0.260031  "
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mean_1             146.718538\n",
       "variance_1           0.590139\n",
       "homogeneity_1        1.607235\n",
       "contrast_1           0.425737\n",
       "dissimilarity_1      0.393991\n",
       "Name: mean, dtype: float64"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out_mean.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "mean_1             0.169601\n",
       "variance_1         0.046337\n",
       "homogeneity_1      0.157245\n",
       "contrast_1         0.098109\n",
       "dissimilarity_1    0.076004\n",
       "Name: std, dtype: float64"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "out_std.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
